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We study the cosmological evolution of domain walls bounded by strings which arise naturally in 
axion models. If we introduce a bias in the potential, walls become metastable and finally disappear. 
We perform two dimensional lattice simulations of domain wall networks and estimate the decay 
rate of domain walls. By using the numerical results, we give a constraint for the bias parameter and 
the Peccei-Quinn scale. We also discuss the possibility to probe axion models by direct detection of 
O |. gravitational waves produced by domain walls. 
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I. INTRODUCTION 
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■ Iii the history of the early universe, various types of phase transitions have presumably occurred as a consequence 
of the spontaneous breaking of symmetries in fundamental physics. During this process, some nontrivial vacuum 
structures, called topological defects, would be formed, and they can affect the evolution of the universe. For ex- 
ample, the formation of cosmic strings has rich phenomenological consequences for cosmology such as the structure 
formations [l[. On the other hand, domain walls, which arise when a discrete symmetry become spontaneously broken, 
, , ■ are considered as cosmologically undesirable objects j2( . 

These topological defects arise in theories beyond the standard model of particle physics. A well-known example is 
£f) . the theory of axions [|| . Axion 0, @ is the pseudo-Nambu-Goldstone boson which arise as a consequence of the Peccei- 
Quinn (PQ) mechanisms @, 0] which is introduced to solve the CP violation problem in quantum chromo dynamics 
(QCD). In this model, one can make the CP violating phase dynamically into zero value by introducing a new 
symmetry called U(1)pq. In the early universe, this symmetry has been spontaneously broken when the temperature 
of the universe has fallen below some energy scale (the PQ scale). This is called the PQ phase transition. At this 
time the networks of global strings called axionic strings are formed. Furthermore, at the time of the QCD phase 
transition, the axion acquires a mass, which causes the spontaneous breaking of discrete Z^ DW subgroup of U(1)pq, 
leading to the formation of domain walls attached to the strings. The integer parameter -/Vdw describes the number 
of degenerate vacua and the number of walls attached to the string. The value of this number depends on the models. 
For example, in the Dine-Fischler-Srednicki-Zhitnitsky (DFSZ) models Nnw — 2iV 9 , where N g is the number 

of quark generations (this can be reduced to Afow = N g , depending on the structure of the Higgs sector [icj|). On 
the other hand, in the Kim-Shifman-Vainshtein-Zakharov (KSVZ) models [U, [Il|], one can get -/Vdw = 1. 

There are two possibilities for the fate of such string- wall networks. One is the case with Now = 1. In this case, 
walls bounded by strings quickly slice themselves into small peaces as discussed in [ID, [l4| • Another possibility is the 
case with A^dw > 1- In this case the networks survive for a long time, which is cosmologically disastrous [15| since 
they overdose the universe and distort the cosmic microwave background (CMB) observed today Q. Therefore, axion 
models with -/Vdw > 1 seem to conflict with the standard cosmology. This is called the axionic domain wall problem. 

There are several ways to avoid this problem. The simplest way is just to consider the model with Now = 1. In this 
model, the walls decay immediately after the QCD phase transition, producing a radiation of barely relativistic axions 
fl6j . This can be realized in the KSVZ models. Also, it is possible to take iVow = 1 in the variant axion models [l7lll8j . 
and its collider implications are recently investigated in 19]. A more intricate solution is to embed the discrete 
subgroup Zn dw of U(1)pq in the center of another continuous group (so called the Lazarides-Shafi mechanism [20l . l2lj |). 
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In this model, the degenerate vacua are connected to each other by another symmetry transformation. However, in this 
kind of model, one have to choose the symmetry group, Higgs representations and U(1)pq charge judiciously, which 
seems to be unlikely to occur (22|. The third possibility is to suppose that inflation has occurred after the PQ phase 
transition to dilute the number density of strings as well as domain walls. In this case, however, a constraint which 
comes from the bound on isocurvature perturbations in CMB observation might be severe [23]. Alternatively, one can 
eliminate domain walls by introducing a tiny Zn dvj breaking term (i.e. a bias) which lifts the vacuum degeneracy 
[HI Ell IS ■ In this case, domain walls collapse due to the pressure force acting between different vacua [25| . 

In this paper we mainly consider the last possibility, biased domain walls, as a solution to the domain wall problem. 
It is natural to expect such a explicit Z./v DW breaking terms in the context of the quantum theory of gravity such as 
string theory: The quantum gravitational effects induce higher-dimensional operators suppressed by powers of the 
Planck mass, which may alter the PQ solution to the CP violating problem I2a-[29ll . It is argued that allowed region 
in the parameter space for this solution is narrow, though it is not ruled out [241 ]. 

The cosmological scenario in which domain walls decay at early time is also interesting from the observational 
point of view. It has been pointed out that the emission of gravitational waves is likely to occur due to the decay 
of domain walls [l6l ]. Indeed, the recent numerical study implies that long-lived unstable domain walls can radiate 
gravitational waves with the amplitude enough to observe [30|. If it is true, it might be possible to probe axion models 
by gravitational wave experiments in the next decades. 

We would like to address these issues more quantitatively, based on the field theory lattice simulations. The 
numerical simulations of axionic domain walls have been performed by several authors. The simulation of decaying 
domain wall with AT DW = 1 was examined in [l6|, |3l|, HH , neglecting the expansion of the universe. Also, the evolution 
of string- wall networks with -ZVdw > 1 models in the expanding background was investigated in [33J. In addition 
to these simulations, we investigate the decay process of the networks due to the explicit Zjv dw breaking term, give 
the numerical confirmation of the decay of the networks with -ZVdw > 1, estimate their lifetime, and constraint the 
parameter space in which the CP violation problem and the domain wall problem are solved simultaneously. 

This paper is organized as follows. In section 2, we describe the model which we consider, and briefly review the 
property of string- wall networks. In section 3, we present the results of the numerical simulations and estimate the 
decay time of the networks. Based on these results, we give observational constraints for the bias parameter and the 
PQ scale in section 4. In section 5, we shortly discuss about a possibility to observe gravitational waves from domain 
wall networks in the future experiments. Then we conclude in section 6. 

Notation: we work in the spatially flat Friedmann-Robertson- Walker (FRW) background in which the metric is 
given by 

ds 2 = dt 2 - a 2 (t)[dx 2 + dy 2 + dz 2 }. 
We denote the cosmic time as t and the conformal time as r, where dr = dt/a(t). 



II. DYNAMICS OF DOMAIN WALL NETWORKS WITH MULTIPLE VACUA 

We consider the model of a complex scalar field (the PQ field) <j> with the Lagrangian density given by 

C = \d^*d»<l>-V{<t>), (1) 

and the potential given by 

V{4>) = ^(0*0 - iff + - cos AW#) + SV, (2) 

4 iV DW 

where 9 is the phase of <fi (i-e. <t> = \(t>\e te ), m corresponds to the mass of the axion, and 77 corresponds to the PQ 
scale. The first term in eq. ^ describes the usual mexican hat potential which causes the spontaneous breaking 
of U(1)pq and produces cosmic strings. Since we are interested in the evolution of domain walls produced after 
the QCD phase transition, we must include the effective potential for the axion field, which is represented as the 
second term in eq. ©. In the absence of the last term in eq. @, the theory has .Zjvdw shift symmetry 9 — » 
9 + 27rfc/A f DW (k = 0, 1, . . . , Now ~~ !)• The spontaneous breaking of this symmetry implies the occurrence of Afow 
domain walls attached to the string. However, we include the additional term in eq. (|2|) which explicitly breaks the 
discrete Zjy DW symmetry [TBI. |24|: 



SV = -£r) 3 ((i>e- iS + h.c). 



(3) 
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Due to this term, domain walls become unstable and eventually decay into the true vacuum at which \S — 9\ is smallest. 
This kind of term may arise from some ultimate theories, but we ignore the detail and simply describe their effect as 
a dimensionless parameter £ (for now, 5 and £ are regarded as free parameters). 

First we consider the case with £ = 0. If we restrict ourselves to low energy configurations at the QCD scale, we 
can put (f) — V el6 an d derive the effective Lagrangian for 8 

£ cS = ^(d^) 2 + ^(l-co S N uw e). (4) 
z "dw 

The equation of motion derived from it in the Minkowski background has the sine- Gordon solution which describes a 
planer wall orthogonal to the z-axis [l[ 

2irk 4 

6{z) = — h — — tan" 1 exp(mz), fc = 0, 1, . . . ,JVnw — !• (5) 

^Vdw A'dw 

From this equation we can estimate the thickness of the wall as 6 W ~ mT 1 . The surface mass density of the wall is 
also estimated as a — 8to?7 2 /./Vq W if we use this exact solution. However, if we include the structure of the neutral 
pion field which varies inside the wall, the surface energy density of the wall becomes a slightly larger value [IH 

a ~ 4.3Um^i/N uw ~ 9mr, 2 /N 2 w , (6) 

where f„ is the pion decay constant and is the mass of the pion. 

By using the exact solution ([5]), we can also estimate the gradient energy stored in the domain wall 

„2 



l lr7jU2 V fd0y 2rj 2 m 2 1 



P^=2\^\ 2 = Y{T Z ) =^coshW 



Note that it is equal to the height of the potential energy 2m 2 r] 2 /N^, w at the center of the wall (z = 0). Therefore, 
if there are planer walls, their gradient energy becomes comparable with their potential energy, as we see in the next 
section. 

It has been argued that the evolution of domain wall networks is characterized by one scale, the Hubble radius, 
and the averaged number of walls per Hubble volume remain the same in the evolution of the universe. Such a 
property is called the scaling solution, and confirmed both numerically [35U38j and analytically (39M4l| for a simple 
model in which the domain walls arise from the spontaneous breaking of Z2 symmetry. Also, the numerical simulation 
performed in [33j indicates that this property is true for networks of iVnw domain walls attached to strings unless 
JVdw = 1. If domain wall networks are in the scaling regime, the energy density of domain walls evolves as 

Pwalls ~ f/t. (8) 

This is equivalent to the fact 

A/Vozt- 1 , (9) 

where A/V is the comoving area density occupied by domain walls. We will check this property in the numerical 
simulation described in the next section. 

If we turn on the bias (£ ^0), walls become unstable and finally disappear. Let us estimate the typical time scale 
for the annihilation of walls. The energy difference between the neighboring vacua introduced by the ^tv dw breaking 
term ([3]) can be estimated as AV ~ 2^ 4 , which acts as a volume pressure py on the wall and accelerates it 
against the false vacuum regions: 

pv~AV~ 47r^ 4 /iVDW- (10) 

On the other hand, the surface tension which straightens the wall up to the horizon scale can be estimated as 

p T - aft. (11) 

The domain walls collapse when these two effects become comparable. From this fact, we can estimate the typical 
time of the decay of domain wall networks as 

tdcc-g/AV- J W 2 , (12) 

where we used eq. ([5]) for a. We will determine the precise value of the numerical coefficient in the formula (1121) from 
the results of the numerical simulations in the next section. 
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III. LATTICE SIMULATIONS 

The dynamical equations which describe the evolution of the PQ field are highly nonlinear. Therefore, we must 
perform numerical simulations in order to obtain any quantitative result. In this section, we describe the results of the 
numerical simulations. First, we present the evolution of unbiased domain walls (i.e. £ = 0) and check the consistency 
of the result with other existing numerical simulations. Next, we investigate the £ dependence of the results and 
estimate the lifetime of domain wall networks. We also comment on the m dependence, which is the dependence on 
the ratio between the axion mass and the PQ scale which we take as a free parameter in the numerical simulations. 
The details for the numerical computations are summarized in appendix [a1 

A. Evolution of the stable domain walls 

We solve the classical field equations derived from the Lagrangian ([T]) [see eqs. (jAll) and (IA2[) ] on the two dimensional 
lattice with 4096 2 points. We used two dimensional simulations since they enable us to investigate the evolution of 
networks with long dynamical range. Memory of the computer limits us to much smaller grid size and shorter 
dynamical ran ge, if we perform three dimensional simulations. In the earlier work about the evolution of the string- 
wall networks [331] . it was confirmed that there are no significant discrepancy between the results of two dimensional 
and three dimensional simulations. 

We choose the value of the parameters as A = 0.1 and 771/77 = 0.1. Note that m/77 represents the ratio between the 
axion mass and the PQ scale. This value for 777/77 is much larger than the actual value of 777/77 by many orders of 
magnitude. It is impossible to use the actual value of 777/77 in the lattice simulations since to -1 represents the width 
of the domain wall, which must be smaller than the box size. Therefore we perform the simulations with 777/77 = 0.1 
and extrapolate the results into the actual value of 777/77. We will mention this point again in section fill B 21 

For stable domain walls, we performed the simulations varying the value of Afow from 1 to 6. The final time of the 
simulation is set to be r/ = 110 in the unit of v^ 1 . Even at this time, the resolution of the wall width and the core 
of stings is well maintained and the Hubble radius is smaller than the simulation box (see appendix I A ip . 

1. Field configurations and time evolution of the area density 

The spatial distribution of the potential energy and the phase of the PQ field are shown in figure [TJ We can see 
that the phase of the PQ field is indeed divided into 7V DW domains and there is the core of strings attached by iV DW 
domain walls at the location where rotates by 2ir. Since these are the results of two dimensional simulations, strings 
exist as a "point" in the two dimensional surface. 

We calculated the time evolution of the comoving area density of domain walls as shown in figure [5J In the case 
with jVdw > 1) the scaling regime in which the area density scales as r _1 [see eq. ([9])] begins around r ~ 10. On the 
other hand, if A^dw = 1, the area density falls off shortly after the beginning of the simulation. These properties are 
similar to that found in the past simulations 33]. 

2. Time evolution of the energy density 

We also calculated the kinetic, gradient, and potential energy densities of the PQ field as shown in figure [31 From 
this figure, we can see that the gradient and potential energy densities evolve as t" 1 cx r~ 2 [cf. eq. ©] from the 
onset of the scaling regime around r ~ 10. However, the kinetic energy decays faster than the gradient and potential 
energies. This is caused by the cosmological expansion that dilutes the energy density of the "radiation" component 
of the PQ field (the oscillation around the minimum of the potential) which redshifts as a~ 4 (t). On the other hand, 
the potential energy and the gradient energy become comparable after the formation of scaling domain walls as we 
have shown in section [TTJ But the gradient energy is slightly larger than the potential energy at the early stage of 
the evolution, especially in the result with large -/Vdw- It might be caused by the fact that if Njyw is large, there are 
many nodes in the networks. In this case, walls are bent in small scale and the gradient energy stored in walls become 
larger than the potential energy. As domain walls straighten up to the horizon scale, the gradient energy becomes 
small and eventually comes up with the potential energy. 

We find some discontinuities in the plots of figure [3j We regard these discontinuities as numerical artifacts rather 
than physical features. The reason is as follows: We note that there is a term proportional to l/\<f>\ in the equation 
of motion of the scalar field [see eqs. (IA1[) and (|A2|) ]. This term becomes singular when the scalar field is on the top 
of the mexican hat potential (\cj>\ = 0) and causes a sudden transition of the field values. This might produce the 
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N DW = 2, x = 056 N DW = 2, x = 056 




FIG. 1: The distributions of the potential energy (left) and the phase of the PQ field 9 (right) with the case JVdw = 2 (top), 
Nnw — 4 (middle), and Now = 6 (bottom) at the time r = 56. The size of these figures is set to be a quarter of the size of 
the simulation box. In the distribution of the energy density, the white region corresponds to the vacuum (V(<f>) ~ 0), the blue 
region corresponds to the domain wall (V(4>) — 2m 2 r/ 2 /N DW ), and the green region corresponds to the string (V(4>) ~ A^ 4 /4), 
but the core of the string is too thin to see in this figure. We take the range of 9 as — tt < 9 < n in the right panel. 

spiky signatures in the plots of the kinetic and gradient energy densities of the scalar field. On the other hand, since 
the potential energy does not contain the term proportional to l/\<j>\ [see eq. (|A3I) ]. there are no discontinuities in the 
plot of the potential energy density. In any case, it seems that this singularity does not affect the evolution of the 
field after the formation of domain walls since the factor l/\cf>\ has a definite value 1/rj. We performed simulations 
in different setups (such as the simulation with the time variable chosen to be cosmic time rather than conformal 
time and the simulations with different spatial resolution) , and found that although the behavior of the kinetic and 
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FIG. 2: The time evolution of the comoving area density of domain walls for various values of JVdw- 



gradient energies in each setup seems to be different, it does not affect the evolution after the formation of scaling 
domain walls. Therefore, we expect that these numerical artifacts hardly alter the final result about the estimation 
of the decay time of domain wall networks. 

In the past numerical simulation performed by Ryden et al. [33j, it was insisted that the ratio between the kinetic 
energy and the total energy of the PQ field tends to be approximately constant and it does not strongly depend on 
the number -/Vdw- However, we can not confirm this fact in our simulations because of the difference in the dynamical 
range: We take the dynamical range as ~ t/ /Tj = 55 in conformal time, but in [33] the dynamical range was taken 
as ~ 200. This is because the simulation performed in [33j relied on the modified equation of motion which keeps 
the resolution of the comoving wall width constant in time to improve the dynamical range of the simulation (called 
PRS algorithm [35]). Now we do not use this formalism since it may erase small-scale structure on the wall [![, which 
might affect the evolution of networks. Therefore, even if the ratio between the kinetic and total energies tends to be 
constant value, it is not seen in our simulation time scale. 



B. Evolution of unstable domain walls and estimation of the decay time 

Next we consider the effect of the term which explicitly breaks the discrete symmetry [i.e. eq. ©]. This effect 
is parameterized by two quantities: 8 and £. In general, 5 can be defined relative to the phase in the quark mass 
matrix and have any value, but in the numerical simulations it only determines the location of the true minimum in 
the potential and does not affect the results of the simulations. Therefore, we take 6 — for a simplicity. 

We note that the lift of the degenerate vacua must be sufficiently small since we assume the circumstance in which 
the discrete symmetry is held approximately. We can understand this requirement more quantitatively as follows. If 
the lift of the degenerate vacua is large, the probability of choosing vacuum at the time of formation of domain walls 
is not uniform between different vacua. For example, assume that there are two vacua (JVdw = 2), and the energy 
density of one vacuum (false vacuum) is greater than that of another vacuum (true vacuum) due to the presence of 
the bias term SV. Then, let us denote the probability of having a scalar field fluctuation at the time of the formation 
of domain walls end up in true vacuum as p t and in false vacuum as pf . The ratio of these two probabilities is given 
by (U 

^=exp(-^Uexp(-i?), (13) 
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where AV^ot — 2m 2 ?7 2 /A r p W is the height of the axion potential, and AVbias — 2£7/ 4 is the difference of the energy 
densities between two vacua. These probabilities are not uniform (i.e. pt = Pj_ = 0.5) if £ 7^ 0. It was shown that 
this non-uniform initial probability distribution also destabilizes domain walls [42I [43| ■ Let us define the parameter e 
which represents the deviation from the uniform probability distribution as p t = 0.5 + e. According to the numerical 
simulations [43J, the time scale of the decay of domain walls Td ec ,prob (in conformal time) due to this effect is given by 

Tdec.prob/Vform — £~ D ^ 2 , (14) 

where Tf orm is the time of the formation of domain walls, and D is the spatial dimension. We must require that this 
time scale should be greater than the simulation time scale (raec.prob > T /)i since we would like to check the effect 
of volume pressure [i.e. the relation given by eq. (|12p] as a decay mechanism of domain walls. In our numerical 
simulations, we choose tj — 110 and Tf orm ~ 10. Then the above requirement leads the condition 



m 2. 



(15) 
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Note that this result is obtained for the case with A^dw = 2. It is not so straightforward to generalize above arguments 
for the case with multiple vacua. However, we guess that the condition for the case with > 2 might be less 

severe than that given by eq. ([IS]) , since if iVow is large, the relative probability of choosing true vacuum becomes 
small and less sensitive to R. Therefore, we take eq. (|15p as a criterion to ignore the effect of non-uniform initial 
probability distribution even for the case with JVdw > 2. It is difficult to satisfy this requirement for large jVdw in 
the numerical simulation since the height of the axion potential is proportional to A r I J w and we have to choose small 
value of £ to make R satisfy the condition (|15p . which is forbidden because of the limitation of the dynamical range 
(domain walls can not decay in the simulation time scale). This restriction forbids us to perform the simulations with 
A^dw = 5 and 6. 



1. Dependence on £ and JVdw 

We performed the simulations for Afow = 2,3,4, and £ satisfying the condition described above and calculated the 
time evolution of the area density of domain walls. The results are shown in figure 2J We confirmed that domain 
walls decay in the time scale which we naively expect as eq. (1121) . 
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FIG. 4: The time evolution of the comoving area density of domain walls for various values of £ with JVdw = 2 (left top), 
JVdw = 3 (right top), and JVdw = 4 (bottom). In the case with JVdw = 4, we can choose only one parameter for £ because of 
the restriction given by eq. (|15[l . 



Now we determine the time scale of the decay of the domain wall networks Td ec based on the results of the numerical 
simulations obtained here. We define Td cc as a conformal time when the value of A/V becomes 1% of that with £ — 0. 
In figure [5] we show the numerical results for Td oc averaging over 5 realizations. We note that Tdcc may depend on 
iVDw as we see i n ec l- d]) ■ Therefore, we plot the dependence on the product iV DW £ in figure [5j We found that the 
results of the simulation with small A/dw£ deviate from the simple scaling, but the results with large A^rjw? can be 
fitted to the naive analytic expectation (|12[) within errors of the numerical results. This deviation in small Adw£ 
might be caused by the finiteness of the simulation box which promotes walls to collapse around the final time of the 
simulation (we observed a similar effect in our previous numerical study [loj]). In order to resolve this effect, we have 
to run simulations with higher spatial resolution (or larger simulation box), which is beyond the scope of this paper. 
Here we just use the results with A^dw£ > 0.001 to fit the data into the relation which we expect from the naive 
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estimation (IT2 



TdccV 



VNvwZ VO.lry 



1/2 



(16) 



We assume that idee is proportional to m as we see in eq. (|12j) and write this dependence explicitly in eq. (|16[) . We 
discuss this point in the next subsection. The best fit value for the numerical coefficient k turns out to be 



k = 2.70 ±0.03. 

From eqs. (|16[) and (j 1 T[) . we get the formula for the decay time of domain walls 



.5 x 



or in cosmic time, 



£dcc — 18 x 



(17) 



(18) 



(19) 



We conclude that the simple scaling of the decay time with the bias parameter is robust at least in the parameter 
range where the finite box effect can be neglected and agrees with simple analytic expectation given by eq. (|12p , except 
the numerical factor of O(10). We will extrapolate this relation for much smaller values of £ in the next section. 
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FIG. 5: The relation between Td cc and Nuw£, obtained from the numerical results with Now = 2 (triangle), A^dw = 3 
(square) and Now = 4 (circle). The dotted line represents the fitting given by eq. I|18[) . which is obtained by using datas with 
Nuwt > 0-001. 



2. Dependence on m 

As we mentioned before, we can not make the value of m/rj arbitrarily small since we invalidate our assumption 
about the approximate discrete symmetry. Furthermore, m/r) can not be so large because the core size of the string S s 
must be smaller than the width of the domain wall 5 W for the formation of stable networks [l6j]: S s /S w ~ m/(r)y/~A) <C 1 
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(for the parameters we used in the simulations, m/{rj\X) ~ 0.32). Therefore, the range of m/rj which we can choose 
in the actual numerical simulation is narrow. 

Figure |6] shows the m/rj dependence of the time evolution of the comoving area density of domain walls. Although 
the range of m/77 is limited, from this figure we can see that the decay time of domain walls (in conformal time) is 
proportional to \frn, as we anticipated in eq. (1121) . Therefore, we assume that this dependence on m is correct for 
other value of m/rj, and use the relation (TT91 for the evaluation of the decay time of domain walls. 
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FIG. 6: The time evolution of the comoving area density of domain walls for various values of m. We set other parameters as 
A = 0.1, iY D w = 2, and £ = 0.0004. 



IV. OBSERVATIONAL CONSTRAINTS 

In the previous section, we treated £ as a free parameter except the condition that it must be sufficiently small. 
However, actual observations can put stringent constraints on this parameter. The purpose of this section is to give 
a quantitative constraints on the parameters by using the numerical result (|T5)) . 



NEDM bound 



The first constraint comes from the fact that the existence of the asymmetric term 8V spoils the PQ solution to 
the CP violation problem in QCD 26-29]. The QCD Lagrangian contains a term 



Co 



32^ 2 



(20) 



where G ap - V is the gluon field strength, G a ^ v = ^e^XcrG aX(J is the dual of it, and 9 is an dimensionless parameter. 
This term violates both P and CP and leads to an neutron electric dipole moment (NEDM) with large amplitude 
which conflict with the experiments. The original PQ idea is to take 9 as a dynamical field with a potential minimized 
at 9 = 0. However, if we introduce the asymmetric term 8V , it shifts the minimum of the potential from 9 = 0. 

Let us define axion field a as <f> = 7je la ^' 1 . Substituting it into eq. ^ and expanding for small a/77, we can derive 
the effective potential for axion field [28[ 
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where m p hy S is the effective mass of the axion 

'phys 

and 6 is the minimum of the potential 



m 2 =m 2 + 2^ 2 cos<5, (22) 



ff = Ai /£\ (23) 
\ T) I m 2 + 2^?7 2 cos o 

The recent upper bound on 8 from the experiments is 9 < 0.7 x 10~ n [H, EH- Assuming £ <C m 2 /rf and sin<5 w 1, 
we have the condition 

- K 2JVrf <07xirl i (24) 

This gives an upper limit on £. 



B. Overclosure bound 



If £ is sufficiently small, string-wall networks continue to exist for a long time. Since they evolve to maintain 
the scaling property as eq. ([5]), their energy density decreases slower than that of matter and radiation, and they 
eventually overdose the universe. If it occurs, they cause the fast expansion and affect the thermal history of the 
universe. Therefore, another constraint comes from the condition that domain walls must not dominate the energy 
density of the universe. 

Let us estimate the time at which domain walls dominate the energy density of the universe. Assuming the 
radiation dominated era, the energy density of the universe and domain walls (at the cosmic time t) are respectively 
p c (t) = (3/8irG)H 2 = 3/327rGi 2 and /9 W aii — a/H^ 1 = cr/2t, where G is the Newton's constant and a is given by 
eq. ([5]). By equating these two, we find that the wall domination occurs at the time 

<wd = idba- (25) 

We must require idee < *DWi from which we obtain the condition 

\ 2 

m 



£>3xl0 3 xA D ^(^-j , (26) 

where we used eq. (fT§|) for i dcc , and M P is the Planck mass. Eq. (121)1) gives a lower bound on £. 

There is a subtlety in this discussion. To make it clear, let us calculate the temperature at which the wall domination 
occurs. By using the relation if^D = (87r 3 ffWD/90A/p)X^, D = 1/(2£wd) 2 , where gwo, Iwd and i?wD are the 
relativistic degrees of freedom, temperature and Hubble parameter at the time i\yD respectively, we found 

/ 1(1 \ i/4 / F \ 1/2 
Twd = 8 x 1(T 2 x ° - MeV. (27) 



VffWD / V!0 12 GeV 

From this equation, we see that the wall domination occurs after the time of big bang nucleosynthesis (BBN) which 
occurs around the temperature Tbbn — lMeV. Furthermore, the upper bound on £ given by eq. (IM1) implies that 
walls must survive until the BBN epoch, as we will see below (in figure [7]) • Therefore, in this scenario string- wall 
networks exist during the BBN epoch, and we must consider whether they spoil the success of the standard BBN 
scenario. 

If the energy density of domain walls becomes sufficiently large during the BBN epoch, it alter the expansion rate 
of the universe, which affects the helium abundance observed today. The energy density of domain walls contributes 
as an extra energy component to the total energy density of the universe during the BBN epoch, which is usually 
parameterized as an effective number of neutrino species N v (see e.g. (45[): 

7T 2 7 

Pextra(iBBN) = -~(JV„ - 3)T BBN . (28) 



Requiring /9 e xtra(*BBN) = /Owaii(^BBN) = vHbbn , where £?bbn is the Hubble parameter at the BBN epoch, we found 

"«- 3 = 5 ■ 8x lJS£= 8 ■ 4xlo M^ , • (29) 
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where c/bbn = 10.75 is the relativistic degree of freedom at the BBN epoch. The recent observations indicate the 
value 3.68-3.80 for N v , but there are various systematic uncertainties [46j. This implies that the extra contribution 
to the energy density at the BBN epoch must be at most N v < 4. From eq. (|29|) we can see that this limit is satisfied 
for all values of F a which do not exceed the cosmological bound [i.e. eq. ([55)) ]. Therefore, domain walls do not have 
a significant effect to the expansion rate during the BBN epoch. 



C. Cold axions from the decay of domain walls 



Domain walls would decay after the BBN epoch and the energy stored in walls is converted mainly into the radiation 
of axion particles and gravitational waves. Let us comment on the effect of the axions produced by domain walls (we 
will discuss about gravitational waves in the next section) . These axions are barely relativistic [l(| [HJ and have a 
Compton length much larger than that of the light elements produced in the BBN epoch. Therefore they can not 
spoil the light element abundance observed today. 

However, since these axions are barely relativistic, they contribute to the energy density of the universe as a cold 
matter component. The energy density of cold axions produced by domain walls must not exceed the energy density 
of the radiations until the time of equality between matter and radiation t eq , which might give a severe constraint for 
this scenario. 

Let us estimate the abundance of cold axions at the time t cq . Suppose that the fraction r of the energy stored 
in the wall become barely relativistic axions at the time of decay of domain walls. Then, the energy density of cold 
axions at the time of decay is given by 

Pardee) = rp wa n(t dcc ), (30) 

where p W aii(^doc) = c/2idcc is the energy density of domain walls at the time t dcc . The number density of cold axions 
at the time t > idee can be estimated as 

Pa(t dec) ( a(idoc) X 3 



uj a V a (*) ) 

where u a = 7w(idcc) is the average energy of the radiated axions, and w(^dcc) is the axion mass at the time tdec- 7 
is the Lorentz factor, and the past numerical study indicates 7 ~ 60 [16]. Using eq. (f3"Tj). we find the fraction of the 
energy density of cold axions at the time t eq as 

„, _ m{t cq )n a (t cq ) _ ra f a(t dcc ) \ 3 
a{ CqJ " Pc(U q ) ~ 2jt cqPc (t cq ) { a(t cq ) J ' {62 > 

where m(t eq ) is the axion mass at the time t oq , and p c (t cq ) is the critical energy density of the universe at the time 
icq- We used the approximation m(t eq ) / 'm(tdec) ~ 1 at the second equality. Substituting eqs. ([5]) and (|19p for a and 
tdec, respectively, and using p c {t cq ) = 2.27 x 10~ 33 (ft M fr 2 ) 4 GeV 4 (see, e.g. [47j]) and a(t dcc ) / a(t cq ) = (V^dccflcq) 1 / 2 , 
where Qm is the ratio of the total matter density today to the critical density, h is the current value of the Hubble 
parameter in units of 100 km sec _1 Mpc _1 and H eq is the Hubble parameter at the time £ cq , we obtain 

Requiring Q a (t cq ) < |, we found the condition 

It gives another lower bound on £. Note that it depends on the value of r, the fraction of the wall energy converted into 
cold axions. If r is as large as O(0. 1), the constraint ([3"4l becomes much severer than that of eq. (f2"6"|). Furthermore, 
it depends on the average Lorentz factor 7 of the radiated axions. The computer simulations performed in 16] found 
7 ~ 7 when \f\r\jm ~ 20. Extrapolating this result into the realistic value \f\r\jm ~ O(10 26 ), they obtained 7 ~ 60. 
However, it is not clear whether this numerical result persists up to the different scale of order 10 26 . If the value of 7 
is actually 0(1) when \/\r)/m ~ 0(1O 26 ), the bound given by eq. (I3"4l becomes much severer, and it may completely 
exclude the scenario which we considering here unless r is considerably suppressed. 
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D. Axion window for a model with -/Vdw > 1 



The constraints (j24|) and ([26l) also depend on the parameter m, which is the axion mass given by 

m ~ 6 x 10" 6 eV ^ 10l ^ GeV ^ , (35) 

where F a = r]/Njyw is the axion decay constant Q. The value of F a has been constrained from various astrophysical 
and cosmological observations. Here we just enumerate axion constraints obtained so far. 

Most astrophysical searches give constraints on the coupling parameters between axions and other particles. Axions 
interact with photons due to the coupling 

■^077 — ^ dF^ F^lu: (^6) 

where F^ is the photon field strength, F^ v = \tp_u\aF Xa is the dual of it, 5,177 = 2 -kF c °77 * s ^ ne axion-photon 
coupling, and a — e 2 /Air is the fine structure constant. The numerical coefficient c a77 is determined by electromagnetic 
and color anomaly and the up/down quark mass ratio z = m u /md (see e.g. Q). Usually it takes a value of 0(1), 
depending on models. Axions also have interactions with fermions 

. dm 



£<W = -i ::2 =- J -aip j "/ 5 ip j , (37) 



where ipj is the fermion field, m,j is its mass, and Cj is a numerical coefficient whose value depends on models. In 
hadronic axion models such as KSVZ model pHll2j. axions do not have interactions with ordinary leptons and quarks. 
On the other hand, in DFSZ model @, Q> axions interact with electrons so that C e = cos 2 /3/3, where tan/3 is the 
ratio of two Higgs vacuum expectation values. 

The axion emission in the sun gives some constraints on the axion-photon coupling. The energy loss by solar axion 
emission requires enhanced nuclear burning and increases solar 8 B neutrino flux. The observation of 8 B neutrino 
flux gives a bound g ail < 5 x 10~ 10 GeV _1 or F a /c ail > 2.3 x 10 6 GeV [48|, |4t|. Also, axion helioscope experiments 
give the bounds on the axion-photon coupling. For example, the recent result of the CERN Axion Solar Telescope 
(CAST) |0| gives a bound g ail < 8.8 x 10" 11 GeV~ 1 , which corresponds to F a /c ail > 1.3 x 10 7 GeV. 

The axion emission from the globular clusters gives another bound on F a . From the helium-burning lifetimes 
of horizontal branch stars, the axion-photon coupling can be constrained as g ail < 0.6 x 10~ 10 GeV _1 (51j, which 
corresponds to F a jc ail > 2 x 10 7 GeV. On the other hand, the delay of helium ignition in low-mass red-giant branch 
stars gives a limit on the axion-electron coupling, a ae e < 0.5 x 10~ 26 or g aee < 3 x 10~ 13 [52?], where a aee = glee/^ 77 , 
g aee = C e m e /F a and m e is the mass of electrons. For DFSZ model, it corresponds to F a /cos 2 f3 > 8 x 10 8 GeV. 

The observation of white-dwarfs also gives the bounds on the axion-electron coupling. The observed luminosity 
of white-dwarfs indicates a aee < 1 x 10~ 26 (53[, which corresponds to F a /cos 2 /3 > 5 x 10 8 GeV for DFSZ model. 
Furthermore, white-dwarfs which show periodic variations in their light curves (known as ZZ-Ceti stars) give the 
more restrictive bound on the coupling a aee . From the observed rate of change of the period of the star G117-B15A, 



it was obtained that a aee < 1.3 x 10 27 or g aee < 1.3 x 10 13 [54]], which corresponds to F a /cos 2 (3 > 1.3 x 10 9 GeV 
for DFSZ model. 

Finally, from the energy loss rate of the supernova (SN) 1987A (55| . one can constrain the axion- nucleon coupling. 
For a small value of the coupling, the mean free path of axions becomes larger than the size of the SN core (so 
called the "free streaming" regime). In this regime, the energy loss rate is proportional to the axion- nucleon coupling 
squared, and one can obtain the limit F a > 4 x 10 8 GeV [49j |. On the other hand, for a large value of the coupling, 
axions are "trapped" inside the SN core. In this regime, by requiring that the axion emission should not have a 
significant effect on the neutrino burst, one can obtain another bound F a < 0(1) x 10 6 GcV 56]. However, in this 
"trapped" regime, it was argued that the strongly-coupled axions with F a < 0(1) x 10 5 GeV would have produced an 
unacceptably large signal at the Kamiokande detector, and hence they are ruled out [571 ] . 

Note that, there is a parameter region around F a ~ 10 6 GeV that cannot be excluded by the SN observation (so 
called the "hadronic axion window") [58[. In particular, in KSVZ model, it is possible that the axion-photon coupling 
vanishes due to the uncertainty of the quark mass ratio m u /md and the bound on g ail described above does not 
significantly constrain the value of F a [58rl60| . In this case, we have to take the hadronic window into account, but 
it might be highly model dependent. In the following, we only take a constraint F a > 4 x 10 8 GeV obtained in the 
"free streaming" regime, since it gives a strong constraint on F a both for KSVZ model and for DFSZ model. We 
take it as a lower bound on F a and neglect other possible model dependences. We note that this lower bound suffers 
from various uncertainties such as the uncertainty of the axion-nucleon couplings and the expression for spin-density 
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structure function which is used to calculate the axion emission rate [49(. These uncertainties may modify the lower 
bound on F a by a factor of 0(1). 

In addition to the astrophysical bounds described above, we must consider the cosmological bound on F a . After 
the QCD phase transition, the axion field behaves like the cold dark matter (CDM) component of the universe due 
to the so called misalignment mechanism [6l| . It is required that the axion CDM abundance must not exceed the 
abundance of CDM observed today. This cosmological requirement gives an upper bound on F a . If we include the 
contribution of axions radiated by string decay, the constraint becomes F a < 3 x 10 11 GeV [62| . 

In summary, we take the astrophysical and cosmological bound on F a as 



4 x 10 8 GcV < F a < 3 x 10 n GeV. 



(38) 



Combining this condition with previous constraints (fM|) . and (|3"4")l . we can draw the allowed region in the £,-F a 
plane as shown in figure [7J From this figure we can see that the allowed region is quite narrow, though it is not 
completely ruled out. We note that there is an unknown factor r. If r is not suppressed, the constraint for the 
model with iVDw > 1 would become much severe, as we see in figure [7] We also note that, as we described above, the 
axion-photon coupling may vanish for a certain value of m u /md in KSVZ model, and the constraint space significantly 
opens up. However, this possibility is not shown in figure [7] 
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FIG. 7: The contours of the various observational constraints in the parameter space of F a and £. The solid line shows the 
bound given by eq. (f2"4")l and the dot-dashed line shows the bound given by eq. (f25|) . for a domain wall number Now = 2 
(red), JVdw = 4 (green) and Now = 6 (blue), respectively. The gray region denoted by "too much NEDM" is the parameter 
space which induces too much CP violation conflicting with the experiment of the NEDM, and the region denoted by "wall 
domination" is the parameter space in which the energy density of domain walls dominates the total energy density of the 
universe. These regions are observationally ruled out. The vertical line shows the usual bound for the axion decay constant, 
given by eq. (|38p . The pink region denoted by "supernovae" is ruled out by the observation of the SN1987A [4j|, and the region 
denoted by "CDM" is excluded by the condition that the axion CDM abundance exceeds the CDM abundance observed today. 
Note that, there are many other astrophysical constraints and various model dependencies (see the text for details). The thin 
solid lines represent the parameters in which Td cc = lMeV and O.lMeV, where Td ec is the temperature at which the domain 
walls decay (we fixed JVdw = 6 for the evaluation of Tdec). The dotted lines show the bound given by eq. (|34[) for r = 0.1 and 
r = 0.01, where we fixed 7 = 60, JVdw = 6 and QmIi 2 — 0.15. The yellow region might be excluded if r is grater than 0.1. 



V. GRAVITATIONAL WAVES FROM AXIONIC DOMAIN WALLS 



If domain walls survive for a long time, they produce gravitational waves. According to the recent numerical 
study [30l | , gravitational waves produced by domain walls have almost flat spectrum from the frequency corresponding 
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to the Hubble radius at idee to the frequency corresponding to the inverse of the wall width, with an amplitude 
determined by a simple dimensional analysis. 

Let us estimate the amplitude and the frequency of gravitational waves produced by axionic domain walls. The 
amplitude of gravitational waves observed today is characterized by a dimensionless quantity f2 gw defined by [45| 

"gw(,roJ = — 7TT— 71 7T> 1^9) 

where p c (to) is the critical energy density of the universe today p gw (to) is the energy density of gravitational waves at 
the present time, and / is their frequency If the gravitational waves are produced at time t„ (assuming the radiation 
dominated era), the amplitude observed today can be written as (63[ 

10 XV3 



n gw (t )h z = 3.60 x io-° I — j n gw (u), (40) 

where g* is the number of radiation degrees of freedom at the time t*, and f2 gw (t*) = Pgw{t*)/pc(t*)- From a 
dimensional analysis, we obtain [30] p gw (£*) ~ Ga 2 , where a is given by eq. (JB]). By using p c (k*) — 3/327rGi 2 and 
t* ~ idec, where idec is given by eq. (fT9"|). we obtain 

nMW * 2 x io- x (l^) 2 (I?)"' 3 (H!|SX) 4 . ,4i, 

The frequency of gravitational waves redshifts as the universe expands. Therefore, the frequency observed today is 
estimated as 

/(io)=^y /3 ^/*, (42) 
\9*J J* 

where /* is the frequency of gravitational waves at the time <*, To =2.725K is the temperature of the universe observed 
today, T* is the temperature at the time t#, and go = 3.36 is the number of radiation degrees of freedom today. The 
numerical result obtained in [30| implies that the frequency of gravitational waves ranges from the Hubble radius at 
the time i* (/* ~ -ff*) to the inverse of the wall width (/» ~ 5~ x ~ m). The frequency corresponding to the Hubble 
radius becomes 

Afo) *3 x 10- x < ? ( T5 i 35 )" J (^W (43) 



where we used the relation 



2 _ 87r 3 g.T, 2 
^* " -WMf ~ 1/(2<dcc) 

to eliminate T* . Similarly, the frequency corresponding to the wall width becomes 

/.<<„> * 5 x 10- x HStf {^ff O" 12 {^f^ <«> 

For example, if we take the parameters Njyw = 4, F a — 10 10 GeV, <?* = 10, and £ = 10~ 58 , which are not ruled out 
by the observations described in the previous section, the amplitude and the frequency of gravitational waves become 
il gw (to)h 2 ~ 5 x 10~ 12 , fh — 2 x 10 _11 Hz, and / ffl ~6x 10 2 Hz. We note that this signal has an amplitude relevant 
to the future gravitational wave experiments. The ground-based experiments are sensitive to the frequency around 
0(10-10 )Hz, which is corresponding to the frequency /„, described above. As a ongoing ground-based experiment, 
LIGO [64| has sensitivities f2 gw /i 2 > 10~ 6 today, and it is planned to upgrade to Advanced LIGO [65| . which would 
have sensitivity C(10 -9 ). The LCGT [66| is also planned with sensitivity similar to that of Advanced LIGO. The 
Einstein Telescope [671 would have more improved sensitivity of O(10~ n ). Also, the space-borne interferometers 
such as LISA [68| and DECIGO [6{| are planned. LISA is sensitive to the frequency around £/(10 _4 -10 _1 )Hz and 
it can reach sensitivity of O(10~ n ). DECIGO is sensitive to the frequency band between LISA and ground-based 
interferometers, and the ultimate sensitivity of DECIGO can reach 0(1O -20 ) [7(3]. Therefore, even if the ground- 
based experiments do not have enough sensitivities, it is possible to detect some signal of gravitational waves in future 
space-borne interferometers. 
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We emphasize that the above argument is based on the numerical result with a simple model in which Z2 symmetry 
is spontaneously broken. Because there is another mass scale 77 (the PQ scale) in addition to the two scales m and 
-ff* , the precise shape of the spectrum of gravitational waves produced by string- wall networks might be different from 
that we obtained numerically in (30j . Therefore, it is necessary to perform three dimensional lattice simulation to 
calculate the spectrum of gravitational waves produced by string-wall networks. Since it is computer intensive task, 
we leave it as a future work. 



VI. CONCLUSIONS 

In this paper, we investigated the evolution of string-wall networks which arise in axion models based on the two 
dimensional lattice simulations. In addition to the usual cosine potential of the axion field, we added the bias term 5V 
which makes domain walls unstable. We confirmed the result of the simulations performed by Ryden et al. [33j : Walls 
with -ZVdw — 1 quickly decay after the beginning of the simulation, and walls with A^dw > 1 evolve maintaing the 
scaling property if we do not include the bias. We also simulated the evolution of unstable domain walls which decay 
due to the existence of the asymmetric term SV in the potential and obtained the time scale of the decay as eq. (|19l) . 
This time scale is consistent with the naive estimation (TT2|) except the numerical coefficient. The observational bounds 
which come from the NEDM and the overclosure of domain walls can constrain the parameter of the bias £ together 
with the axion decay constant F a as eqs. (f2~4"|) and ([26"]) , and it is found that the allowed region is rather narrow. In 
addition to these constraints, we must require that the abundance of cold axions produced by the decay of domain 
walls should not exceed the CDM abundance at the time of equality between matter and radiation. This gives another 
bound (|34p . and it depends on the fraction r of the wall energy converted into cold axions. This constraint might be 
severer than that obtained from the condition of the overclosure of domain walls unless r is considerably suppressed. 
However, if such string-wall networks have existed, gravitational waves produced by domain walls would be detectable 
in future interferometers. 

This issue provides a new observational method for axion physics in the sense that we might be able to probe 
axion models via gravitational wave astronomy. If a signal of gravitational waves which we described in section [V] is 
detected in the future experiments, it might be possible to measure the parameters such as a bias £ and the axion 
decay constant F a from the spectrum of gravitational waves. Even if we can not detect any signal, we may constrain 
the window of the allowed region described in section HVDI or completely rule out the model with A^dw > 1< 

For now we estimated the amplitude and the frequency of gravitational waves based on the simple dimensional 
analysis, but the precise shape of the spectrum of gravitational waves is not determined. Especially, the spectrum 
of gravitational waves produced by string-wall networks which arise when Z^ DW symmetry is spontaneously broken 
might be different from that with the simple model in which Z2 symmetry is spontaneously broken. Also, the value of 
r, which is the fraction of the wall energy converted into cold axions when the decay of walls, is undetermined. The 
constraint for the axion models with A^dw > 1 might become much severer if the value of r is not so small. These 
undetermined factors will be settled by direct numerical simulations, which we will tackle in the future work. 
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Appendix A: Notes on the numerical computations 

In this appendix we describe the setup for our numerical computations. 
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Formulation 



We define two real scalar fields 4>\ and <f>2 such that <f> = 0i + ife- In the numerical studies, we normalize the 
dimensionful quantities in the unit of r\. For example, — > (f>/f], t — > tr], etc. With this normalization, the equations 
of motion for two real scalar fields derived from the Lagrangian ([!} with the FRW background are 



V 2 , . , ,. 2 , , 2 . _ _ m 2 



i + 3i?0i - 7277T01 = -A0i(0| + 2 I - l) + 2£cos<5+ — — — sin 6> sin N DW 9, (Al) 



a 2 (i)^ ™ x ^ ' * JV DW |0| 

and the potential energy can be written in the form 

A iM- ].)- -!-- 

V DW 



^2 -7277a ^2 = - A0 2 (0f + 01 - 1) + 2g cos (5 - m _ JM cos sin A DW 0, (A2) 



^(0) = t(I0| 2 - l) 2 + T3-(l - cos A DW #) - 2|0|£cos(0 - 5) + 2£, (A3) 



where we add a constant for V(0) in order to make V(0) = at the minimum of the potential. 
In the numerical simulations, we use the conformal time r, and normalize by the scale factor, 

d _ 1 d ^_ 4> 
dt a dr a ' 

4> + 3tf = -I U" - —i] , (A5) 



(A4) 



where a prime represents a derivative with respect to r. Note that, in the radiation dominated universe, a" /a = in 
the last equation. Then, from eqs. fjAip , (|A2|) and (|A5[) , we obtain 

4 2 

— — — — CJ 77? 

0' 1 '-V 2 1 = -A0i(|0| 2 -a 2 ) + 2a 3 £cosc5 + ^ sinflsin A DW 6», (A6) 

A DW |0| 



4 2 

^'-V 2 2 = -A0 2 (|0| 2 - a 2 ) + 2a 3 gcos<5 - - m , -, cos6>sin N uw d. (A7) 

A DW |0| 



We choose the initial time of the simulation so that t% = 1 and a(tj) = 1. This corresponds to the condition tj = 2 in 
conformal time. 

The simulations are performed in the comoving box with size b (in the unit of ?y _1 ). The lattice spacing is Sx — b/N 
where N is the number of grid points (here we take N = 4096). The physical scale of the Hubble radius H^ 1 = 2t, 
the width of the wall 5 W = to -1 , and the core size of the string S s = 1/vA divided by the physical lattice spacing 
fephys = a(t)Sx are respectively 

H- 1 Nt <L, N /T4\ , 5, N 



(?)■ - itz-M)- < A8 » 



fephys b (fcphys bm \T ) fephys 

In our simulation we take b — 230, m — 0.1 and A = 0.1. Then, at the end of the simulation t/ = 110, these 
ratios become H~ x /Sx p h ys — 1959 < N, 5 W /Sx p h ys ~ 3.2, and 6 s /6x p hy S — 1-02. Therefore, these length scales are 
marginally resolvable even at the final time of the simulation. 

We put the periodic boundary condition in the configuration of the scalar fields. We solve the time evolution of the 
fields by using the fourth order Runge-Kutta method. 



2. Initial conditions 



The initial conditions are similar with that used in [71] . We treat 0i and 02 as two independent real scalar fields 
so that each of them has quantum fluctuation at the initial time with correlation function in the momentum space 
given by 

<&(k)&(k')) = ^(27r) 3 <J< 3 >(k + k'), (A9) 

<&(k)&(k')) = |(27r) 3 ^(k + k'). (i = l,2) (A10) 
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Since the effective squared masses of the fields quickly become negative at the initial time (i.e. on the top of the 
mexican hat potential), we used the massless fluctuations as the initial conditions, replacing the factor \/k 2 + to 2 
with k in the above formulae. We also put the momentum cutoff fc cu t above which all fluctuations are set to zero in 
order to eliminate the unphysical noise which comes from high frequency modes in the field distributions. Here we 
set fc cu t = 1. 

Note that, although we describe the results of the computation on the "two dimensional lattice" in the body of 
paper, the dimensionality of delta functions in eqs. (|A9|) and (|A10[) is 3. We performed numerical simulations on the 
two dimensional lattice, which means that we choose the number of lattice point in z-axis to be N = 1 while we take 
N = 4096 in x- and y- axes. This choice corresponds to a three dimensional field theory with plane symmetry along 
z-axis. 

We generate initial conditions in momentum space as Gaussian random amplitudes satisfying eqs. (| A9|) and (|A10I) . 
then Fourier transform them into the configuration space to give the spatial distributions of the fields. It is possible 
to start the simulation with more realistic condition, for example, with the thermal initial conditions. Thermal initial 
conditions are relevant if we are interested in the formation of global strings, since the formation time and initial spatial 
distribution of strings are affected by the initial amplitude given at high temperature. After the string formation, 
the field distribution just before the QCD phase transition is determined by the vacuum configuration around global 
string networks. Then, the formation of axionic domain walls occurs if the expansion rate of the universe becomes 
less than the mass of axion. Therefore, thermal initial conditions are not relevant in the QCD era. In any case, it 
is difficult to reproduce the whole process described above in the numerical simulations because of the limitation of 
the dynamical range. However, since it seems that the scaling property is not so much affected by the initial field 
configurations and we are interested in the evolution of the fields after the formation of the string-wall networks, we 
expect that the results were qualitatively unchanged if we used the different initial conditions. 



We used the similar algorithm for calculation of the area density of domain walls with that used in [30j. We sum 
up the quantity S± which takes either 1 or for each of the links (the neighboring grid points that differ by one) with 
the weighting factor as 



where 6, x , etc. is a derivative of #(x) with respect to x, and C is a coefficient chosen so that A/V = 1 is satisfied 
when all the links have the value 5± = 1. If -/Vdw = 1) we put S± — 1 at the location where fa has different signs and 
</>i < at the two ends of the link (i.e. at the location on which 9 — n). If Now > 1, we divide the codomain of 9 
(2ir) into Now domains, and put S± = 1 at the location where 9 belongs to different domains at the two ends of the 
link. Otherwise we put 5± — 0. 
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